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Abstract 

We call a system bouncing ball billiard if it consists of a particle that is subject to a 
constant vertical force and bounces inelastically on a one-dimensional vibrating pe- 
riodically corrugated floor. Here we choose circular scatterers that are very shallow, 
hence this billiard is a deterministic diffusive version of the well-known bouncing 
ball problem on a flat vibrating plate. Computer simulations show that the diffusion 
coefficient of this system is a highly irregular function of the vibration frequency ex- 
hibiting pronounced maxima whenever there are resonances between the vibration 
frequency and the average time of flight of a particle. In addition, there exist irreg- 
ularities on finer scales that are due to higher-order dynamical correlations pointing 
towards a fractal structure of this curve. We analyze the diffusive dynamics by clas- 
sifying the attracting sets and by working out a simple random walk approximation 
for diffusion, which is systematically refined by using a Green-Kubo formula. 

Key words: fractal diffusion coefficient, bouncing ball, granular material, phase 
locking 

PACS: 47.20.Ky, 46.40.Ff, 47.53.+n, 45.70.-n 



1 Introduction 



An important challenge in the field of nonequilibrium statistical mechanics is 
to achieve a more profound understanding of transport processes by starting 
from the complete microscopic, deterministic, typically chaotic equations of 
motion of a many-particle system [1,2,3]. However, up to now such a detailed 
analysis concerning the microscopic origin of nonequilibrium transport could 
only be performed for certain classes of toy models. Very popular among these 
models are systems that are low-dimensional, spatially periodic and consist of 
a gas of moving point particles that do not interact with each other but only 
with fixed scatterers. 
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For these chaotic dynamical systems it was found that, typically, the respec- 
tive transport coefficients are highly irregular functions of the control pa- 
rameter. These irregularities were exactly calculated for diffusion in a simple 
abstract piecewise linear map [4,5,6] and were shown to be related to topo- 
logical instabilities under parameter variation. The same properties are exhib- 
ited by diffusive Hamiltonian particle billiards such as the periodic Lorentz 
gas [7,8], the flower-shaped billiard [9], and for a particle moving on a one- 
dimensional periodically corrugated floor [10]. Similar irregularities occur in 
nonhyperbolic maps with anomalous diffusion [11] in simple models exhibiting 
(electric) currents [12,13] as well as in the parameter dependence of chemical 
reaction rates in reaction-diffusion processes [14]. In most of the cases men- 
tioned above, the deterministic transport coefficients clearly showed fractal 
structures [4,6,8,10,11,12,14]. Seemingly analogous transport anomalies were 
measured in experiments [15,16] or have been reported for models that ap- 
pear to be rather close to experiments [17,18]. However, still an experimental 
verification of the fractal nature of irregular transport coefficients remains an 
open question. 

In this paper we wish to contribute to this problem by introducing and inves- 
tigating a dynamical system that, we believe, is rather close to specific exper- 
iments. For this purpose we consider a generalized version of the interesting 
model described in Ref. [10], where a particle subject to a constant vertical 
acceleration makes jumps on a periodically corrugated floor. We generalize 
this model by including some friction at the collisions and compensate this 
energy loss by periodic oscillations of the corrugated floor. This generalized 
model is an example of what we call a bouncing ball billiard. It is designed to 
be very close to the experiments on diffusion in granular material performed in 
Refs. [19,20,21], thus we hope that our theoretical predicitions can be verified 
by respective experiments. However, we emphasize the fact that, in contrast 
to most experiments on granular material,^ here we study a granular gas 
that consists of statistical ensembles of a one-particle system only. One should 
therefore be very careful in concluding anything from the physics we discuss 
here for the case of highly interacting many-particle systems. 

The corrugated floor of our model is formed by circular scatterers that are de- 
liberately very shallow. Hence, another important limiting case of our model 
is the famous bouncing ball problem, where an inelastic particle bounces verti- 
cally on a flat vibrating surface. This problem has both been studied experi- 
mentally [22,23,24,25,26] as well as theoretically [27,28,29,30]. It is well-known 
that the nonlinear dynamics of the bouncing ball is very complex exhibiting 
one or more coexisting attractors depending on the frequency of the floor [31], 
that is, the system displays both ergodic and nonergodic dynamics. The main 



1 An exception is the experiment performed in Ref. [19] where the dynamics is 
exclusively due to the collisions of a single particle with a corrugated surface. 
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theme of this paper will be to show that dynamical correlations in terms of 
phase locking, as associated to these different attractors, have a huge impact 
on diffusion in the bouncing ball billiard and that these regimes lead to strong 
irregularities in the parameter-dependent diffusion coefficient on large and 
small scales. 

The paper is composed of six sections: In Section 2 we outline some important 
features of the bouncing ball problem. This knowledge provides a roadmap in 
order to understand the diffusive behavior in the spatially extended system. 
In Section 3 we give the full equations of motion of the bouncing ball billiard 
and briefly explain the specific choice of the control parameters. In Section 4 
we present results for the diffusion coefficient and for the average energy of 
the system as functions of the frequency of the vibrating floor. In Section 5 we 
analyze the irregularities of the parameter-dependent diffusion coefficient in 
full detail by relating them to different dynamical regimes, as represented by 
certain structures in the corresponding attracting sets. In Section 6 we work 
out a simple random walk approximation somewhat explaining features of the 
diffusion coefficient on the coarsest scales. This approach is then refined by 
including higher-order dynamical correlations based on a Green-Kubo formula 
for diffusion thus providing detailed evidence for the existence of irregularities 
on finer and finer scales. Section 7 contains a summary of the results, an outline 
of further links to previous works, and an outlook concerning further studies 
in this direction. In particular, we are trying to encourage an experimental 
realization of our bouncing ball billiard. 



2 The bouncing ball problem: Arnold tongues 

A standard problem of chaotic dynamics that was extensively studied is the 
one of a ball subject to an external field that bounces inelastically on a flat 
vibrating surface [22,23,24,25,26,27,28,29,30,32]. At first view this model ap- 
pears to be rather simple, however, its equations of motion are in fact highly 
nonlinear and generally do not allow any exact solution. Hence, a good part 
of the investigations deals with simplified versions of this model such as the 
so-called high bounce approximation leading to the dissipative standard map 
[26,31]. Here we restrict to the exact model only as nicely analyzed in Refs. 
[26,30] and recall some important features that will be needed to understand 
the forthcoming results. 

Let us assume that the flat surface exhibits a sinusoidal motion y = —A sin(wt), 
where A and oo are the amplitude respectively the frequency of the vibration. 
Between the bounces on the surface the ball moves in a gravitational field with 
acceleration g. If the mass of the ground plate is assumed to be much larger 
than the mass of the ball, the latter becomes a trivial quantity in the equa- 
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tions of motion that, for convenience, we set equal to one. The bouncing ball 
problem then has four control parameters: A, u, the gravitational acceleration 
g and the vertical restitution coefficient a. After a proper transformation of 
the equations of motion A, u and g can be grouped into T = Au 2 /g. A par- 
ticular complexity of the bouncing ball system is related to the fact that for 
given values of the parameters V and a the system may posess more than one 
attractor, each with a different basin of attraction. This is drastically exempli- 
fied in a study of the high-bounce approximation revealing the coexistence of 
an arbitrarily large number of attractors if the damping is small enough [31]. 




Fig. 1. The three major Arnold tongues of the bouncing ball problem in parameter 
space. The adjacent odd and even lines represent the boundaries of the resonances 
between the period of the driving T v ib r and the time of flight between two bounces 
Tfjight, see Eq. (1). From bottom to top it is Tfli g ht/T v ib r = 1,2,3. a and T are 
dimensionless quantities. 

We will now focus on the Arnold tongues of the bouncing ball defining a series 
of stable resonances. The largest ones were analyzed approximately in Refs. 
[23,25] and were lateron calculated exactly in Ref. [30]. An account of the 
hierarchy of smaller ones in the limit of zero friction was approximately given 
in Ref. [29]. A resonances of order kj\ means that the period of the particle's 
flight between two collisions Tfl ig ht is k times larger than the period of the 
floor's motion, T flight = kT vihr = k 2ti/uj. The region in parameter space (a, T) 
where a resonance of order k exists can then be calculated via linear stability 
analysis of the equations of motion to [30] 
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These different major tongues are shown in Fig. 1. Note that, somewhat in con- 
trast to simple Arnold tongues, here not all the inital conditions in phase space 
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may exhibit phase locking. That is, typically there are in addition nonresonant 
trajectories that will "stick" to the surface interrupted by sequences in which 
they are relaunched exhibiting jumps of smaller and smaller amplitude. This 
particularly complicated dynamics was denoted as "self-reanimating chaos" 
[25] or as "chattering" [30]. In Ref. [30], a careful analysis indeed led to the 
result "that generic trajectories starting under experimental conditions ter- 
minate in a region of chattering", which seemed to be at variance with the 
experimental observation of a Feigenbaum-like bifurcation scenario and a re- 
spective transition to chaos [22,23,24,25,26]. 



3 The bouncing ball billiard 



3. 1 Equations of motion and control parameters of the model 
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Fig. 2. Illustration of the bouncing ball billiard studied here: A moving point particle 
is colliding inelastically with circular scatterers that are periodically distributed on 
a line. Parallel to y there acts an external field with constant acceleration g. The 
floor of scatterers oscillates with an amplitude A and a frequency to. 

The bouncing ball billiard that we study in this paper, with the floor formed 
by circular scatterers, is depicted in Fig. 2. 

The equations of motion of this system are defined as follows: The particle 
performs a free flight in the graviational field g \\ y between two collisions. 
Correspondingly, its coordinates (x~ +1 ,y~ +1 ) and velocities (v~ n+1 ,v~ n+1 ) at 
time t n+ i immediately before the n + 1th collision and its coordinates y+) 
and velocities (v£ n ,Vy n ) at time t n immediately after the nth collision are 
related by the equations 



X n+1 ~ X n + V xn(tn+1 tn) 

Vn+l =Vn+ Vynttn+l ~ *„) ~ ~ fn)7 2 > 



11 == 1] 

x,n+l u xn 

11 == 11~^~ 

y,n+l yn 



g(t n +i ~ tn) ■ 



(2) 
(3) 
(4) 
(5) 
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At the collisions the change of the velocities is given by 



vf n - v ciLn = a (v ci±n - v ±n ) (6) 

V{ n ~ Uci|| n = /5(ujf n -Vd\\n) , (7) 

where v ci is the velocity of the corrugated floor. Here we distinguish between 
the different velocity components relative to the normal vector at the surface 
of the scatterers, where the scatterers are represented by the arcs of the cir- 
cles forming the floor. v±, v\\ and v C i±, tyy indicate the normal, respectively 
tangential components of the particle's, respectively the floor's velocity with 
respect to the surface at the scattering point. Correspondingly, we introduce 
two different restitution coefficients a and (5 that are perpendicular, respec- 
tively tangential to the normal. 

As in case of the standard bouncing ball problem we assume that the floor 
oscillates sinusoidally, y ci = —A sin (cut), where A and u are the amplitude 
respectively the frequency of the vibration, see Fig.2. Hence, the phase space 
is defined by the variables (t,x,y,v x ,v y ), where the time variable needs to 
be included because we have a driven system described by nonautonomous 
differential equations. If we record only the collision events with the floor, as 
it is usually done in case of the bouncing ball, the dimension of the phase space 
can be reduced accordingly. For this purpose we use the position of the collision 
at the circumference of a scatterer defined by the clockwise arclength at the 
impact s and the horizontal and vertical velocities just after a collison, z = 
(s,v£,Vy). In such a Poincare section the time-continuous flow is equivalent 
to the time-discrete map 



Zn+l = /(z n ) 

t n+ i =t n + r(z n ) (8) 

ln+1 = In + a ( z n) ■ 

called a suspended flow in Refs. [1,10,33]. In the second equation r(z) is the 
time of flight between two successive collisions, and in the third equation 
a(z n ) counts the number of scatterers, or tiles, the particle has jumped over 
during a flight, l n being the tile where the particle has actually started at 
the nth collision, hence l n is an integer. If the billiard is simple enough, these 
equations can be derived analytically, such as in case of the periodic Lorentz 
gas [1,33] and in case of the completely elastic non- vibrating floor [10]. This 
enormously simplifies the computer simulations since the mapping has merely 
to be iterated numerically. However, in case of the bouncing ball billiard due 
to the vibrations we were not able to derive the exact form of this mapping. 
Hence, the equations of motion were solved via discretizing the parabolic flight 
between collisions in time and trying to accurately determine the point of the 
collision with the surface by reducing the time lag of the iterations to 10~ 6 s 
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at collision events. 



3.2 Choosing the parameters of the billiard 

We choose the parameters of the billiard in Fig. 2 according to a possible 
experimental setup [21]. The extension of each scatterer we determine to d = 
2mm. In our study we fix the value of a and increase the value of T. This can be 
done by tuning the frequency / = u/2tt while keeping the amplitude constant 
at A = 0.1mm. The gravitational acceleration is taken to be g = 9.8m/s 2 . 
The value of the normal restitution coefficient is chosen to a = 0.5 such that 
the major resonances as depicted in Fig. 1 do not overlap with each other. 
Since our model has the geometry of a dispersing billiard, the arcs induce a 
dynamical instability that somewhat counteracts the stability in presence of 
a resonance. More importantly, the arcs make the system diffusive by feeding 
energy that is vertically pumped into the system also into the respective x- 
component of the velocity. 

The main purpose of our model is to study the impact of dynamical correla- 
tions on the deterministic diffusion coefficient, hence we would like to preserve 
the mechanism related to the resonances of the bouncing ball problem as much 
as possible.0 Consequently, the radius of the scatterers should be very large. 
Here we choose the radius R of the arcs and the corresponding curvature to 
K = 1/R = 0.04mm -1 . Note that, because of the spatial periodicity, this 
billiard is of the form of a two-particle problem with periodic boundary con- 
ditions. Since the radii of the scatterer and of the particle are simply additive, 
a periodic surface with respective radii could be realized experimentally by 
choosing a suitably large radius of the moving particle (which, of course, may 
have other side effects on which we do not elaborate here). 

At parameter values where no resonance is possible the particle may "stick" to 
the surface for some time, that is, it will land at a certain position and, because 
of the friction, it will be relaunched at the next period of the oscillation. 
However, in case of fully elastic tangential scattering with (3=1, because of the 
spatial extension long trajectories may occur at which a particle "slides" along 

2 We briefly remark that, by starting from the curved surface of Ref. [10] and 
making the originally Hamiltonian dynamics dissipative, first simulations seemed 
to indicate that the irregular structure of the diffusion coefficient was rather im- 
mediately destroyed. That may be related to the fact that certain periodic orbits 
associated to small islands of stability, which indirectly determined this structure, 
are straightforwardly eliminated by inclusion of dissipation. Furthermore, the oscil- 
lations of the plate profoundly disturbed the original Hamiltonian dynamics. The 
precise nature of the transition between the Hamiltonian case and the dissipative 
case may be an interesting problem for further studies. 
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the surface in one direction even if its vertical velocity relative to the surface 
of the scatterer is zero. In order to avoid these rather unphysical horizontal 
slides, we set the value of (5 smaller than one. Numerically it turned out that 
a horizontal restitution of f3 — 0.99 is already sufficient to eliminate all slides. 
We remark that this value is still considerably larger than the ones found in 
experiments on inelastic particle collisions [34], however, with the choice of 
this value we avoid the vanishing of the diffusion coefficient at certain values 
of the frequency. 

Instead of using the dimensionless parameter Y we make our following presen- 
tations in terms of the frequency /, since with a variation in steps of 0.2Hz it 
affects T = AAtt f 2 /g only in the third digit after the comma, consequently 
the frequency provides a much better scale. This translates the position of the 
Arnold tongues in the bouncing ball problem at a = 0.5, as shown in Fig. 1, as 
follows: 1/1-resonance: / G [50.98, 61.56]Hz; 2/1-resonance: [72.10, 76. 71]Hz; 
3/1-resonance: [88.30, 90. 95]Hz. The position of the tongues is also indicated 
in Fig. 3. In the following, the diffusion coefficient is presented in units of 
mm 2 /s and the energy in kgmm 2 /s 2 . 



4 The frequency-dependent diffusion coefficient of the bouncing 
ball billiard 

If we choose the set of parameters as explained in Sec. 3.2, we expect that 
the dynamics of the spatially extended system is to a large extent somewhat 
determined by the one of the bouncing ball problem. However, as was outlined 
in Sec. 2, for the bouncing ball problem time and ensemble averages often 
strongly depend on the inital conditions, i.e., the system is nonergodic. In the 
bouncing ball billiard this deficiency is eliminated by the defocusing character 
of the scatterers, i.e., the dynamics is getting ergodic, except for some short 
frequency intervals at small frequencies around the 1/1-resonance. That is, the 
dynamics typically evolves on only one or sometimes on two attractors. The 
details related to these phase space properties will be worked out in Sec. 5. 
Here we focus on the coarse structure of the diffusion coefficient as related to 
regions of resonances. 

For any frequency located in the interval of / G [49.77, 95. 0]Hz the system 
was found to evolve into a nonequilibrium steady state. This was checked in 
terms of the average (kinetic, total) energy of the moving particle, which was 
always approaching a constant value for large times. The diffusion coefficient 
in the horizontal direction was then obtained from the Einstein formula 




(9) 
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where the brackets denote an ensemble average over moving particles. The 
diffusion coefficient was converging to a specific value for large enough times. 
The size of the ensemble of points consisted of 10 3 initial conditions that 
uniformly filled the phase space region of < s < 2.0004mm, \v£\ < 40mm/s, 
and < v+ < 40mm/s.Q 

The dependence of the diffusion coefficient on the frequency is depicted in 
Fig. 3. This figure clearly demonstrates that whenever the value of the fre- 
quency reaches the different major Arnold tongues as discussed in Sec. 2, the 
diffusion coefficient has considerably larger values than elsewhere and conse- 
quently exhibits local maxima. 

Above the frequency regime corresponding to the 1/1 Arnold tongue, that is, 
for / > 70Hz, the diffusion coefficient becomes very small but is still different 
from zero in that it has a value around one or two. This phenomenon lasts 
until / ~ 80Hz. Note that for a smaller tangential restitution coefficient this 
region would collapse to D — 0mm 2 /s. Fig.3 shows that at the frequencies 
associated to the 1/1 and to the 3/1 resonances the diffusion coefficient is 
about two orders of magnitude larger than the ones at / € [70, 80]Hz. 

Irrespective of the specific frequency value, after a transient time the average 
of the total energy converges to a constant value indicating the existence of a 
nonequilibrium steady state. These values are presented in Fig. 4 as a function 
of the frequency. This curve resembles very much the parameter-dependent 
diffusion coefficient. However, as may be expected intuitively the value of the 
diffusion coefficient is even more closely associated with the average kinetic 
energy related to the horizontal velocity component t^/2, see Fig. 5. Still, in 
detail the structure of the diffusion coefficient cannot be trivially understood 
on the basis of the frequency dependence of the (kinetic) energy only, or vice 
versa; see Sec. 6 for further details. It is aslo interesting to note here that there 
is no equipartitioning of energy between the x and the y coordinate. Fig. 4 
furthermore includes a lower bound for the energy of a moving particle, which 
presents the case when it sticks to the surface and behaves like a simple har- 
monic oscillator. This approximation indicates an increase of the total energy, 



For a considerable number of frequency values the results for the diffusion coef- 
ficient have been cross-checked by sampling the initial conditions for the ensemble 
from one long trajectory. That is, we take as initial conditions coordinates from 
this trajectory after the dynamics has evolved a random fraction of a maximal time 
that is considerably larger than the time of typical transients. Within errorbars, the 
same results have been obtained for the diffusion coefficient. However, this cross- 
check could not easily be performed for frequency intervals marked by (b)(i) in Fig. 
3 where the dynamics is nonergodic. In these cases the ensemble was still the one 
described in Sec. 4, apart from the two points indicated in Fig. 3, and the results 
may thus depend on the choice of initial conditions if other attracting sets were not 
sampled by the specific set of initial conditions. 
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Fig. 3. The diffusion coefficient of the bouncing ball billiard as a function of the 
vibration frequency of the corrugated floor. The frequency range of the non-diffusive 
bouncing ball resonances, see Fig.l at a = 0.5, are shown on top of the frame. They 
suggest a strong impact of the resonances on the diffusion coefficient. However, 
the fine structure of this curve appears to be more complex than merely explained 
by the resonances. Different dynamical regimes, as thoroughly discussed in Sec. 5, 
are indicated in this figure in form of horizontal arrows. The two crosses (x) at the 
frequencies 74.0Hz and 74.5Hz reveal the existence of the second resonance but were 
only obtained by sampling a larger set of initial conditions. The standard deviation 
errors are smaller than the magnitude of the symbols. 

and of the corresponding horizontal kinetic energy, on a large scale which, 
however, appears to be largely obscured by the superimposed oscillations of 
the energy. 

A further interesting feature of these three figures is that above / ~ 80Hz 
particularly the diffusion coefficient, and to some extent also the average en- 
ergies, seem to experience a pronounced average increase on larger scales. This 
property somewhat reminds of the phase transition-like behavior discussed in 
Ref. [19] where the average horizontal velocity of a granular multilayer of par- 
ticles in an asymmtric periodic potential shows a sudden growth from zero to 
a nonzero value. 

Finally, we would like to draw the attention to a particularly interesting detail 
of these curves: Around / ~ 60Hz in the first resonance there appears to be an 
almost perfect linear decrease of the average total energy, and correspondingly 
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Fig. 4. The average total energy as a function of the vibration frequency. The dotted 
line shows the possible minimal energy a particle can have, which is the energy of 
a harmonic oscillator sticking to the surface, A 2 u> 2 /2. 
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Fig. 5. The average horizontal kinetic energy, that is, Eki n ,x = v%./2, as a function of 
the vibration frequency. One can see that the structure of the kinetic energy closely 
follows the one of the diffusion coefficient in Fig. 3. 

of the associated kinetic energy, whereas just in the same region the diffusion 
coefficient exhibits some rather regular oscillations on fine scales. One may 
argue that here the variation of the energy, as a function of the frequency, 
provides a linear "ramp" that essentially reduces the bouncing ball billiard 
to the system of Ref. [10], where the diffusion coefficient was studied as a 
function of the energy as a control parameter. Indeed, the diffusion coefficient 
of Ref. [10] too exhibited some very regular, apparently fractal oscillations 
under variation of the energy, which qualitatively very much resemble the 
ones that appear in the respective small region of Fig. 3. We will get back to 
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a more refined analysis of the fine structure of the whole curve in Sect. 6. 



5 The different dynamical regimes of the bouncing ball billiard 

Depending on the driving parameter T, or respectively on the frequency /, the 
following dynamical regimes can be distinguished in the bouncing ball billiard: 

(a) When T < 1 the maximal acceleration of the floor is smaller than the 
gravitational force, Auj 2 < g, consequently once the particle lands on the 
surface it will never leave it again. In this case the particle is eventually located 
in one of the wedges that the arcs of two adjacent circles form with each other. 
Hence, after a short transient period the particle will be in phase with the 
surface and its relative velocity to the floor will be zero. For our parameters 
this happens up to the frequency / = 49.77Hz. For frequency values just 
above /o a particle can be relaunched from the floor, however, because of the 
inelasticity of the collisions and due to the height of the arcs typically the 
particle remains trapped in the wedges of the billiard. Numerically we find 
that around / = 57Hz particles start to leave a wedge for the first time, as 
shown in Fig. 3, consequently we consider this value as the onset of diffusion, 
despite the drastic increase that takes place right afterwards around / = 58Hz. 

(b) The most interesting dynamical regime appears to be associated with the 
1/1-resonance of the bouncing ball. In case of a flat vibrating plate there are 
two attractors at these frequencies, one showing this resonant behavior and 
another one with particles that stick to the surface. In case of the bouncing 
ball billiard this dynamical regime exhibits a complex scenario under varia- 
tion of the frequency, where due to the convexity of the scatterers also ergodic 
motion becomes possible. This suggests to distinguish between two different 
types of dynamics: (i) There are frequency regions where the system has two 
coexisting attractors, one representing the 1/1-resonance, and a second one 
showing an intermittent-like behavior consisting of long periods of stick and 
slip, or "creepy" motion. This nonergodic regime is located at the frequency 
intervals / G [57.0, 57.8]Hz, / G [60.4, 61.8]Hz and / G [63.4, 65.6]Hz, see Fig. 
3. (ii) The other type of behavior is represented by frequencies at which there is 
only one attractor with all inital conditions leading to the same phase locking. 
These frequency intervals are at / G [58.0, 60.2]Hz and at / G [62.0, 63.2]Hz. 
It is remarkable that whenever particles are in resonance with the plate the 
motion becomes regular only in the vertical direction. In this direction we 
then have a 1/1 phase locked dynamics similar to the respective resonance of 
the bouning ball. However, in the horizontal direction we still find a highly 
irregular behavior as quantified by the diffusion coefficient. In order to illus- 
trate this dynamics we look at a projection of the phase space at the collisions 
of a particle with the floor. Because of the resonance the vertical velocity 
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has approximately the same value at the collision such that we almost have 
a Poincare section of the dynamics. Thus we only consider the position s of 
the colliding particle on the circumference of a scatterer and the horizontal 
velocity v x right after a collision, see Fig. 6. 

60 I 1 1 1 1 1 1 1 1 1 1 
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02 . 
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0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 
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Fig. 6. Projection of the phase space in the 1/1 resonance where the vertical dynam- 
ics is approximately in phase with the vibration frequency, which is here / = 58Hz. 
s marks the clockwise position of the colliding particle at the circumference of a 
scatterer and i>+ denotes the horizontal velocity right after a collision. 

In this figure one can detect traces of a separatrix structure around the un- 
stable fixed point at s — l,v£ — corresponding to the two cases where a 
particle can leave a wedge, i.e., the trajectory goes over s = 1, or it cannot, 
i.e., the particle does not approach the peak of an obsctacle or simply turns 
back in the neigbourhood of s = 1. In addition, there are interesting parts of 
the phase space outside the separatrix in the regions of s < 1, i>+ > 0, respec- 
tively s > 1, v£ < 0. Here particles go "uphill" by generating a fan-shaped 
structure in phase space showing that specific paths are chosen according to 
which particles are allowed to leave a wedge, whereas other possibilities are 
apparently forbidden. Yet there is no evidence that particles creep along the 
surface of a scatterer in form of long sequences of very tiny hops, cp. to Fig. 7, 
however, one may speculate that the fan-shaped structure presents a precur- 
sor of such creeps. We furthermore remark that for this figure the dynamics 
is ergodic and that the average energy of a particle saturates after several 
hundred collisions such that we get a good convergence at this frequency. 

(c) Between the frequencies / = 65.8Hz and / = 81.0Hz the dynamics is char- 
acterized by long creeps. Consequently, leaving the wedges takes a long time 
and suppresses diffusion. In the neighbourhood of the frequency 74Hz the 2/1- 
resonance could be detected, however, this resonance could only be revealed 
by choosing a set of initial conditions that is considerably more extended 
over the phase space than at the other frequencies, with initial conditions of 
v + > lOOm/s. The results indicated by the two specially marked points in Fig. 
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3 at the frequencies 74Hz and 74.5Hz show that the impact of this resonance 
on the diffusion coefficient is, according to the respective measure of such or- 
bits in the nonergodic phase space, rather small. Note that, apart from this 
special case, we had no indication for other non-ergodic behavior in regime 



We furthermore note the existence of another small peak around 77Hz that is 
still close to the region of this resonance. Curiously, the following simple rea- 
soning precisely identifies this frequency as a special point of the dynamics: 
Let us assume that a particle behaves similar to a harmonic oscillator stick- 
ing to the surface. Then, as argued above, its average vertical kinetic energy 
is equal to E y = A 2 uj 2 /2. However, for diffusion this energy should be large 
enough such that a particle can pass an obstacle, which is here of the height 
h = 0.02mm. Consequently, if the particle aquires a potential energy that is 
equal to E pot = g(A + h) it can pass any scatterer even at the highest position 
of the vibrating plate. Of course, in this simple reasoning we have disregarded 
any impact of the two restitution coefficients as well as the fact that, addition- 
ally, some horizontal kinetic energy is necessary in order to perform a diffusion 
process. In any case, a respective calculation yields / = 77.23Hz precisely cor- 
responding to the position of this small peak. One may furthermore speculate 
that this argument is somewhat related to the onset of the global increase of 
the diffusion coefficient and of the energy as discussed in Sec. 4. 

(d) For the next, in this study the last dynamical regime with / e [81.1, 95.0]Hz 
Figs. 4 and 5 demonstrate a noticeable increase of the average total and hori- 
zontal kinetic energies that appears to be reflected in a respective increase of 
the diffusion coefficient on a coarse scale. Now the particle spends considerably 
more time for moving around than sticking to the surface. Indeed, we find that 
the dynamics is ergodic such that, in contrast to region (b), here both types 
of resonant and creepy motion are intimately linked to each other, instead 
of appearing in coexisting attractors, or at different frequencies. This type of 
dynamics is represented by the attractor of Fig. 7, where the coordinates s 
and v x are shown at the freqency / = 85.2Hz. 

In spite of the fact that particles creep from time to time along the surface ,Q 
as we can clearly see in the attractor of Fig. 7, the average energy of a particle 
launched from randomly chosen initial conditions on this attractor saturates. 
Consequenty, a system of noninteracting particles reaches a nonequilibrium 
steady state in this regime, too. 

At the frequency interval of the 3/1 resonance one can now find a quasi- locked 
orbit, see Fig. 8, that still induces a strong maximum in the diffusion coefficient. 
At first view this orbit indeed appears to resemble a period 3 orbit, however, 

4 For the probability and the respective fraction of time of sticking to the surface 
see [30,32]. 
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Fig. 7. Projection of the phase space at the frequency of / = 85.2Hz representing a 
dynamical regime where particles "creep" along the surface. These creeps are visible 
in form of the almost vertical lines indicating that a particle performs long sequences 
of correlated jumps nearly at the same position with smaller and smaller horizontal 
velocities. The coordinates are the same as in Fig. 6. 
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Fig. 8. A sequence of a particle's trajectory at / = 88.4Hz (left) and the correspond- 
ing attractor in the plane (y,v+) (right). These figures show the deformation of the 
previous 3/1-resonance into a much more complicated quasi-periodic orbit. 
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closer inspection of the y coordinate reveals that it may rather have bifurcated 
into portions that resemble an orbit of period 6. The right half of Fig.8 depicts 
these six regions that are visited one after the other. However, an even closer 
look reveals that in fact this is a period 11 orbit: the trajectory first follows the 
positions 1 to 6, but at any second sequence it skips the position 6 and goes 
instead directly from position 5 to the position 1. Even more, sometimes the 
particle misses the last landing at position 4 and gets stuck at the surface. It 
is then relaunched at the next period, with the next bounce after the relaunch 
at position 5, and the periodic orbit starts again. 
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In summary, for all frequencies there appear to exist essentially four differ- 
ent types of dynamical regimes: if the diffusion coefficient is very large, the 
dynamics typically exhibits some type of periodicity that is responsible for 
an enhancement of diffusion. Either this dynamics is strictly associated to 
an original bouncing ball resonance, as in case of the regimes (b)(ii), or the 
resonance still appears in a deformed way as in case of the largest peak in re- 
gion (d). If the diffusion coefficient is getting smaller there is an intermediate 
regime in which two different attractors with resonances and creepy motion 
coexist, and the dynamics is non-ergodic, see regimes (b)(i) and the small 
peak around / = 74Hz. However, it is also possible that the motion in this 
intermediate regime is intermittent-like and ergodic, by alternating between 
localized creeps and long flights. Finally, if the diffusion coefficient is getting 
closer to zero, there is typically only sticky motion as discussed for region 
(c). Fig. 6 exemplifies the regime of large diffusion coefficients, whereas Fig. 7 
shows portions of typical creepy dynamics. Note that, since we do not know 
about the position of any higher-order Arnold tongues in the phase diagram 
Fig. 1, we cannot possibly associate smaller peaks in the diffusion coefficient to 
such higher-order tongues, but we cannot exclude that this is possible. In any 
case, in the next section we outline a simple approach that helps to explain 
the origin of the structure of the diffusion coefficient particularly on smaller 
scales. 



6 The fine scale of the coefficient: correlated random walk approx- 
imations 



As we explained in the previous section, the huge peaks in the diffusion coeffi- 
cient are associated with the major resonances of the bouncing ball problem. 
However, first of all we have not yet clarified what large scale functional de- 
pendence one would expect for the diffusion coefficient based on stochastic 
theory. In addition, one observes a lot of irregularities on finer scales in Fig. 3 
whose origin needs to be commented upon. In this section we work out some 
simple random walk arguments as well as systematic refinements of it, which 
were called correlated random walk approximations [8], in order to achieve a 
more detailed understanding of Fig. 3. 

We first start with a simple approximation that is based on the picture of 
diffusion as a random walk on the line. Let us assume that the wedges of the 
billiard act as "traps" for a bouncing particle in which it stays for an average 
time r. Let us furthermore assume that, after this time, the particle escapes 
just to the neighbouring trap to the left or to the right. Disregarding any 
correlations between subsequent jumps the random walk diffusion coefficient 
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is given by 



DrM) = 



2r(/) ' 



(10) 



where d is the distance between two wedges and r is the escape time of a 
particle out of a wedge. For simple billiards such as Lorentz gases and related 
models this formula can be worked out exactly [9,8,35]. But here the dynamics 
is more complicated, hence we first obtain r from computer simulations by 
restricting ourselves to frequency regions in which the dynamics is ergodic. 
The diffusion coefficient Eq. (10) resulting from this numerical input is shown 
in Fig. 9. 



q 20 




80 

/(Hz) 

Fig. 9. Numerical and analytical evaluations of the random walk approximation Eq. 
(10) for the diffusion coefficient with frequencies / > 65.6Hz. The bold line shows 
again the diffusion coefficient of Fig. 3, the dotted monotonously increasing curve 
represents the analytical formula Eq. (13), and the dashed irregular curve close to 
zero depicts Eq. (10) with the escape time computed numerically. The analytical 
approximation indicates some increase of the diffusion coefficient on coarse scales, 
whereas the numerical evaluation yields an irregular structure that is qualitatively 
very similar to the one of the precise diffusion coefficient. 

The irregularities resulting from D rw are qualitatively very close to the ones 
of the precise diffusion coefficient. This is easily understood by approximating 
the average escape time to 



d 

<v x > 



d 

\r2Ex 



(11) 



hence the frequency dependence of the escape time must be very closely related 
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to the frequency dependence of the horizontal kinetic energy E x as shown in 
Fig. 5, whose functional form was already observed to be close to the one of 
the diffusion coefficient. 

Eq. (11) furthermore provides a convenient starting point for a simple analyt- 
ical approximation of the random walk diffusion coefficient Eq. (10). In order 
to find an analytical estimate for E x we start from the energy balance 



where E is the average total energy and E y gives the average vertical kinetic en- 
ergy. E pot = gy is the average potential energy of the particle. We approximate 
the average height y with the amplitude of the vibrations, y ~ A, and the aver- 
age total energy with the energy of the harmonic oscillator, E ~ A 2 uo 2 /2. The 
conjecture that the total energy of a moving particle is roughly proportional to 
u 2 is also supported by a simplified version of a vibrating plate [36]. Finally, 
from computer simulations we approximate the relation between horizontal 
and vertical kinetic energy to E y ~ 19E X reflecting the fact that only a small 
percentage of the energy vertically pumped into the system is transferred into 
the horizontal direction. This appears to be due to the very shallow geometry 
of the corrugated floor and to the action of the tangential friction. Taking all 
these expressions into account we get from Eq. (12) 2E X ~ A 2 u 2 /20 — gA/lt). 
Inserting this expression into Eq. (11) and combining the resulting formula for 
r with Eq. (10) we get a simple analytical result for the random walk diffusion 
coefficient which is 



This formula is depicted in Fig. 9. After a discontinuous, phase transition-like 
onset of diffusion it indicates some systematic but slow increase of the diffusion 
coefficient on a coarse scale. In the region of / e [80, 90] there is an analogy 
of this increase to the one obtained from the numerical evaluation of Eq. (10). 

We now systematically refine the simple approximation Eq. (10) in order to 
further explore the origin of the irregular structure of the diffusion coefficient 
on fine scales. For this purpose we apply the scheme proposed in Ref . [8] , which 
suggests to approximate the precise diffusion coefficient by adding higher-order 
terms to Eq. (10) according to a Green-Kubo formula for diffusion. 

Let us still assume that a particle moves with a frequency of jumps 1/r from 
traps to traps situated on a one-dimensional lattice at a distance d. Starting 
from the Einstein relation Eq. (9) one can prove that in the limit of infinitely 



E = E x + E y + E pot , 



(12) 




(13) 
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Fig. 10. Magnification of the frequency interval / £ [84.6, 86.4]Hz in Fig. 3. The 
bold line on top shows the respective irregular diffusion coefficient on a fine scale, 
the lowest dotted line depicts the corresponding random walk approximation. The 
other lines represent a series of higher-order approximations refining this random 
walk and converging to this curve. These approximations D n (f) are obtained from 
the truncated Green-Kubo formula Eq. (15) for n = 0, 2, 4 and 7. The standard 
deviation errors are smaller than the magnitude of the symbols. 

many steps Eq. (9) is equal to the sum 

1 oo 

D(f) = -Y,Ck<h(x )-h(x k )>, (14) 

T k=0 



which is a version of the Green-Kubo formula for an ensemble of particles 
diffusing in a billiard [8]. Here h(xk) defines jumps of a particle at position x k 
at the kth time step in terms of respective lattice vectors, which here take the 
values ±d. The first coefficient is cq — 1/2, and for k > it is c k — 1. 

We now truncate this series after the nth time step yielding the systematic 
sequence of diffusion coefficient approximations 

D n (f) = f + - p(s lS2 ...)h-h( Sl s 2 ...). (15) 

^ T T s 1 s 2 ...s n 

Here Sis 2 ■ ■ ■ denotes symbol sequences suitably identifying a particle's tra- 
jectory passing through the traps [8]. The averages in Eq. (14) are now given 
in terms of the conditional probabilities p(sis 2 ■ ■ ■) that characterize all these 
possible trajectories. The first term in Eq. (15), that is, D (f), is again the 
random walk approximation Eq. (10). 
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Let us now sketch of how to evaluate Eq. (15) to higher order. Let the first 
step be in the positive (right) direction with h(xo) = d. The next jump is then 
defined by h(xi) = d si , where si can be either / (forward) when the particle 
jumps to the right (in the same direction as the first one) or b (backward) 
when the particle jumps to the left, with df — d and db = —d. If the dynam- 
ics is deterministic and highly correlated, a particle starting from some tile 
may return to it with a probability that is different from the backscattering 
probability of a random walk, which is 1/2. Hence, by feeding in the precise 
value for this probability as, e.g., computed from simulations one can suitably 
correct the random walk diffusion coefficient D (f) in first order. 

In detail, the first approximation D\{f) at time step 2r reads 



depending on the backscattering probability pb, respectively the forwardscat- 
tering probability pf, only. In our case the value of pf is considerably larger 
than the value of Pb, which is due to the fact that the height of the arcs of 
the circles between two traps is very shallow, thus yielding a considerable 
contribution to D\. In a similar way the second correction can be evaluated 
to 



The probabilities associated to different symbol sequences were all computed 
numerically, and for the frequency interval of / e [84.6, 86.4]Hz the resulting 
approximations of increasingly higher order are shown in Fig. 10. Evidently, 
the random walk diffusion coefficient D (f) is quantitatively much smaller 
than the numerically obtained value of the diffusion coefficient. However, by 
adding more and more terms representing repeated forward- and backward- 
scattering this series of approximations converges to D(f). Note that in or- 
der to obtain a good approximation one has to go to relatively high order. 
That way, this figure exemplifies the impact of strong dynamical correlations 
on the diffusion coefficient. Furthermore, this convergence is obviously not 
monotonous in the frequency, with the higher order approximations revealing 
irregularities on fine scales. This behavior is well-known from related systems 
exhibiting deterministic diffusion [4,5,6,7,8,9,10,11] and typically indicates the 
existence of a diffusion coefficient that is highly irregular as a function of the 
parameter, or possibly even fractal. The origin of these irregularities on finer 
scale may be attributed to the geometry of the scatterers, which make the 
system topologically instable and induce a severe pruning of trajectories that 
is reflected in the diffusion coefficient. In other words, under parameter varia- 
tion certain orbits may be forbidden, but the nature of forbidden orbits may 
change in a very intricate way [37]. 



D l = D + 2D (p f - Pb ) =D + 2D (1 - 2p b ) 



(16) 



D 2 = D l + 2D (p f f - p fb + Pbf ~ Pbb) ■ 



(17) 



20 



In summary, we analyzed the frequency dependent diffusion coefficient on 
different scales: On a coarse scale and for large enough frequencies the diffusion 
coefficient appears to increase in qualitative agreement with a simple analytical 
random walk approximation. However, this behavior is largely suppressed by 
a series of huge peaks that we attribute to the impact of resonances, which 
can be traced back to phase locking in the bouncing ball problem. Finally, 
successive evaluations of the Green-Kubo formula of diffusion give evidence 
for the existence of further irregularities on finer and finer scales that indicate 
higher-order dynamical correlations, and one may speculate that this curve is 
even of a fractal nature. 



7 Conclusions 

The purpose of this work was to analyze deterministic diffusion in a novel 
class of systems that we denoted as bouncing ball billiards. The geometry of 
this class of models is very similar to the one of some well-known chaotic 
Hamiltonian single-particle billiards such as Lorentz gases [1,2] in consisting 
of a periodic array of scatterers together with a point particle that collides 
with these obstacles. A particularly interesting case was studied in Ref. [10] 
composed of a periodically corrugated one-dimensional floor in which a par- 
ticle experiences a constant vertical acceleration. Here evidence was given for 
the existence of a highly irregular, possibly fractal diffusion coefficient, as it 
was already found to exist in some very related, but somewhat simpler chaotic 
maps [4,5,6]. Our goal was to move this mechanical system closer to recent 
experiments on granular material, where similar geometries have been used. 
Hence, we made the collisions inelastic by introducing normal and tangential 
restitution coefficients. In this sense our model is somewhat related to the in- 
elastic Lorentz gas of Ref. [38] . However, in order to compensate for this loss 
of energy at the collisions the surface was moving according to an oscillating 
driving force, just as it was done in respective experiments. Finally, we com- 
posed the periodic surface of circular scatterers that are very shallow. This 
way, our system is very close to the well-known problem of a ball bouncing 
inelastically on a flat surface, which is highly nonlinear and exhibits strong 
dynamical correlations in terms of phase locking. Hence, this setup enabled us 
to study features of chaotic diffusion in a specific granular one-particle system. 

Our main numerical result is the existence of a highly irregular, possibly frac- 
tal diffusion coefficient as a function of the driving frequency. We performed 
a detailed analysis in order to understand the complicated structure of this 
curve and found that there are many different dynamical regimes depending on 
the driving frequency. The most pronounced effects are due to some previous 
regions of phase locking of the bouncing ball problem that still manifest them- 
selves in form of huge peaks in the diffusion coefficient. Further features of our 
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dynamics are parameter regions in which the particle's dynamics may "stick 
to" and "creep along" the surface. Although mostly associated to a very small 
value of the diffusion coefficient, this dynamics may also come together with 
a quasi-periodic orbit again exhibiting a huge peak at the position of a former 
bouncing ball resonance. Finally, creeps and resonances may also coexist in 
form of some nonergodic dynamics evolving on two different attractors. We 
furthermore detected signs of a phase transition-like behavior at higher fre- 
quencies that, to some extent, seemed to match to a simple analytical random 
walk approximation. Refined approximations in form of correlated random 
walks, based on the systematical numerical evaluation of a Green-Kubo for- 
mula of diffusion, yielded evidence for the existence of further irregularities 
on finer and finer scales pointing towards a possible fractal structure of the 
frequency-dependent diffusion coefficient. 

Due to the fact that both the bouncing ball and diffusion of granular par- 
ticles on periodic surfaces are systems that have been profoundly studied in 
experiments [19,21,22,23,24,25,26], we hope that this paper stimulates further 
experimental work in order to verify this irregular frequency dependence of the 
diffusion coefficient. Of course we would not expect to reveal any fractal curve 
in a real experiment. However, we do believe that at least the largest peak is 
rather stable against random perturbations, as appears to be corroborated by 
the studies in Refs. [39,40]. Consequently, phenomena like the 1/ 1-resonance 
may show up in form of a local maximum of the diffusion coefficient as a func- 
tion of the frequency. Again, we emphasize at this point that without further 
studies regarding the impact of perturbations on this billiard we cannot con- 
clude to which extent dynamical correlations such as the ones reported in this 
paper may survive in case of a gas of interacting many particles. However, 
we note that in Refs. [20,41] the 1/1-phase locking regime of the bouncing 
ball was made responsible for the experimentally observed [20,42] phase tran- 
sition between a quiescent amorphous and a gaseous state in a monolayer of 
a granular gas. 

Of course, more theoretical work needs to be done in order to predict more 
details of possible experiments. There already exists an experimental setup 
measuring diffusion of single granular particles in one dimension [19], which is 
very close to our present model, however, the respective channels are presently 
asymmetric. On the other hand, computer simulations of a two-dimensional 
vibrating plate with circular scatterers would be close to the experimental 
device of Ref. [21]. Furthermore, of course there exist more realistic models 
for inelastic collisions than by using two constant restitution coefficients [34]. 
For example, one may wish to include the possibility of transfer of rotational 
energy of the moving particle at the collision. 

One may also want to investigate the impact of a tilt on the bouncing ball 
billiard making it into a ratchet-like device. On the basis of our studies and 
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related to the findings of Ref. [12], one may suspect that the current may ex- 
hibit again a highly irregular, possibly fractal structure as a function of con- 
trol parameters reflecting the intrinsic properties of highly correlated chaotic 
transport. First signs of such irregularities indeed appear to be present in the 
numerical and experimental data of Ref. [19]. It may also be elucidating to 
compute the spectrum of Lyapunov exponents and possibly other dynamical 
systems quantities, as well as time-dependent correlation functions and prob- 
ability distributions of position and momentum variables for such types of 
bouncing ball billiards^ and to check their characteristics under parameter 
variation. Although simple reasoning suggests that the bouncing ball billiard 
is chaotic whenever there is a resonance, in form of providing at least one posi- 
tive Lyapunov exponent, the features of the dynamical instability in parameter 
regions of stickiness is far from clear. 

Another crucial point is to learn more about the dynamics of the standard 
bouncing ball problem with respect to the existence of higher-order resonances. 
As was shown in Ref. [29] , for an approximate bouncing ball map in the limit 
of zero friction coefficient the complete set of tongues forms a very complicated 
fractal Devil's staircase-like structure in the parameter space. Unfortunately, 
for the general case up to now only the major Arnold tongues have been cal- 
culated, and to find more tongues, to adequately construct more details of 
the respective phase diagram, and to check in more detail for the existence of 
bifurcations in the inelastically bouncing ball problem appear to be open ques- 
tions. However, in our case such information provided a road map in order to 
predict, and to understand, the irregular structure of the frequency-dependent 
diffusion coefficient, hence a respective analysis would be fundamental for fur- 
ther research along these lines. 

Finally, we wish to remark that similar studies regarding the impact of phase 
locking on diffusive properties have been performed in Refs. [17,44]. However, 
in both cases the systems under investigation were very different from the 
granular model studied here: Ref. [17] outlined the impact of phase locking on 
the magnetoresistance in antidot-lattices under constant electric and magnetic 
fields and equipped with a bulk friction coefficient [17], whereas in Ref. [44] 
the authors analyzed the impact of phase locking on the diffusion coefficient 
for an overdamped stochastic Langevin equation with a ratchet-like asymmet- 
ric periodic potential under the influence of a periodic tilt force. The same 
Langevin equation, however, with a symmetric potential was used in Ref. [45] 
revealing again an enhancement of diffusion due to the matching of an average 
waiting time of a particle with the driving frequency. This nice paper as well as 



5 Note that periodic fine structures of probability distributions that are analogous 
to the ones reported in Ref. [43] were also detected in the chaotic maps of Ref. [5], 
hence we would conjecture that the bouncing ball billiard shows the same type of 
distributions representing basic features of deterministic dynamics. 
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the very related studies of Ref. [46] point to a rather rich literature regarding 
the possible impact of stochastic resonance-like phenomena on the diffusion 
coefficient. Of course, the bouncing ball billiard too exhibits more important 
time scales, of which the average escape time is an example, hence we would 
indeed conjecture that there exist further resonances with the driving period 
in this model that may lead to additional modulations of the diffusion coef- 
ficient. To find such parameter regions further adds to the list of interesting 
open problems related to the bouncing ball billiard. 
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